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We examine the stability of a class of solitons, obtained from a generalization of the Boussinesq 
equation, which have been proposed to be relevant for pulse propagation in biomembranes and 
nerves. These solitons are found to be stable with respect to small amplitude fluctuations. They 
emerge naturally from non-solitonic initial excitations and are robust in the presence of dissipation. 
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INTRODUCTION 

The action potential in nerves is a propagating volt- 
age pulse across the axonal membrane with an ampli- 
tude of about 100 mV. In 1952, A. L. Hodgkin and A. 
F. Huxley proposed a theory for the nerve pulse which 
has since become the textbook model Q. Their pic- 
ture is based on the equilibration of ion gradients across 
the nerve membrane through specific ion-conducting pro- 
teins (called ion channels) which leads to transient volt- 
age changes. Hodgkin-Huxley theory thus relies on dis- 
sipative processes and is intrinsically not isentropic. It 
is rather based on Kirchhoff circuits involing capacitors 
(the nerve membrane), resistors (the ion channels) and 
electrical currents introduced by the ion fluxes. 

We have recently proposed an alternative model for 
the nerve pulses based on the propagation of a local- 
ized density pulse (soliton) in the axon membrane |2j. 
This model has several advantages over the Hodgkin- 
Huxley model. It explains the reversible temperature 
and heat changes observed in connection with the nerve 
pulse. (Such reversible changes are not consistent with 
the Hodgkin-Huxley theory but rather suggest that an 
isentropic process is responsible for the action potential 
[J El-) It further predicts the correct pulse propa- 
gation velocities in myelinated nerves. These velocities 
are closely related to the lateral sound velocities in the 
nerve membrane. One essential feature of our model is 
the presence of empirically known lipid phase transitions 
slightly below physiological temperatures. The closer the 
phase transition is to physiological temperatures, the eas- 
ier it is to excite the nerve pulse. Our model therefore 
immediately explains another interesting feature of nerve 
excitation, i.e., that the nerve pulse can be induced by a 
sudden cooling of the nerve,and that it can be inhibited 
by a temperature increase During compression, the 
appearance of a voltage pulse is merely a consequence of 
the piezo-electric nature of the nerve membrane, which 
is partially charged and asymmetric. 

Another advantage of a soliton-based description of 
pulse propagation in nerves lies in its predictive power. 
Given measured values of the compression modulus as a 
function of lateral density and frequency, soliton prop- 



erties (including its shape and its energy) can be deter- 
mined uniquely as a function of soliton velocity. Given a 
measured soliton velocity, the theory contains no freely 
adjustable parameters and has the virtue of being falsi- 
fiable. 

In the possibility of soliton propagation was ex- 
plored and compared to observations in real nerves. In 
the present paper we study some intrinsic features of 
these solitons, in particular the stability of such pulses 
in the presence of noise and dissipation. Such investi- 
gations are necessary to demonstrate that such pulses 
could propagate under realistic physiological conditions 
over the length scales of nerves (as much as several me- 
ters) even in the presence of friction and lateral hetero- 
geneities. In the following section, we will state the model 
more precisely and derive the analytic form of its soli- 
tonic solutions. We will then turn to a description of 
the numerical methods used here. We use these methods 
to probe (i) the stability of solitons with respect to "in- 
finitesimal" perturbations (i.e., lattice noise), (ii) the way 
in which solitons are produced by localized non-solitonic 
initial excitations of the system, and (iii) the behavior of 
solitons in the presence of dissipation. We will demon- 
strate that the solitons of Ref. [2| are remarkably robust 
with respect to all of these perturbations. 



ANALYTIC CONSIDERATIONS 

Thermodynamic measurements of the lipids of biolog- 
ical membranes reveal a number of interesting features 
of potential relevance for understanding the nature of 
pulses in biomembranes and nerves. In particular, such 
systems display an order-disorder transition at tempera- 
tures somewhat below that of biological interest from a 
low temperature "solid-ordered" phase to a high temper- 
ature "liquid-disordered" phase in which both the lateral 
order and chain order of the lipid molecules is lost 
The proximity of this phase transition to temperatures 
of biological interest has striking effects on the compres- 
sion modulus and, hence, on the sound velocity [9! Ho|. 
For densities some 10% above the equilibrium density, 
the low-frequency sound velocity is reduced by roughly 
a factor of 3 from the velocity of cq = 176.6 m/s found 
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at equilibrium. The sound velocity then rises sharply, 
returning to the value Co at a density roughly 20% above 
the equilibrium density. Measurements at high frequen- 
cies (i.e., 5 MHz) reveal a much smaller dip in the lateral 
compression modulus and a sound velocity that is always 
materially larger than that at low frequencies and thus 
indicate the presence of significant dispersion 0, 0] . 

In Ref. 0] , these features were exploited to suggest that 
the propagation of sound in these lipid mixtures can be 
described by the equation 
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- p A is the difference between the lateral 



mass density of the membrane and its empirical equilib- 
rium value of p A = 4.035 x 10~ 3 g/m 2 , and the low fre- 
quency sound velocity is cq — 176.6 m/s. The coefficients 
p and q were fitted to measured values of the sound ve- 
locity as a function of density. Although high frequency 
sound velocity measurements indicate that the dispersive 
coefficient, h, must be positive, neither the magnitude of 
h nor the specific form of this term have been verified ex- 
perimentally. In practice, the only role of h is to establish 
the linear size of solitons, and it can thus be chosen, e.g., 
so that the width of the soliton is comparable to that 
known for nerve pulses. Here, we choose to work with 
the dimcnsionless variables u, x and t defined as 
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With this choice of variables, eq. QJ assumes the form 
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The qualitative features of the empirical compression 
modulus require that B\ < and B 2 > 0. In the 
numerical work described below, we will adopt the pa- 
rameter values B\ = —16.6 and B 2 — 79.5 found in 
Ref. 0. Eq.© can be recognized as a generalization of 
the Boussinesq equation, and it is known to have expo- 
nentially localized "solitonic" solutions which propagate 
without distortion for a finite range of sub-sonic veloci- 
ties. We now determine the analytic form of these soli- 
tons. 

Since we seek solutions which propagate without dis- 
tortion, we regard u as a function of £ = x — fit and 
rewrite Eq. J2J as 



B 2 — = — 



(5) 



We can integrate this equation twice with the assumption 
that u vanishes at spatially infinity to obtain 
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It is clear from this equation that exponentially localized 
solutions are possible if 1 < 1. Multiplication by du/d^ 
and a final integration leave us with the result 
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It is clear that u is symmetric about its maximum value. 
The solution will grow from until it reaches a maximum 
value at which du/d^ = 0. Equation @ indicates that 
this is possible only if 
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For the parameters B\ = —16.6 and B 2 = 79.5 adopted 
in 0|, we find /3q ~ 0.649851. We will use these parameter 
values in the remainder of this paper. 
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FIG. 1: Soliton profiles for velocities /3 = /3 +4 x 10 , 0.65, 
0.734761, 0.85, and 0.95. The maximum height diminishes as 
a function of j3. The width of the soliton diverges for both 
/3 — + f3o and /3 — > 1 and has a minimum at (5 « 0.734761, 
which corresponds to the dashed curve. 

We thus expect localized solutions for /3q < |/3| < 1. 
When this condition is met, the right side of Eq. (J7J will 
have two real roots, u — a± with 




(9) 



It is readily verified that the desired solitonic solutions 
of Eq. © have the analytic form 
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These solutions are shown in Figure ^ f° r a selection of 
soliton velocities. 

As expected, Eq. © can be obtained from a suitable 
energy density. We thus seek an energy density, £, such 
that Eq. will result from variation of the correspond- 
ing Lagrangian density. To this end, it is useful to intro- 
duce the dimensionless displacement, s(x, t), defined as 
u = ds/dx. The energy density can then be written as 
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The two terms in Eq. (|llf) represent the kinetic and po- 
tential energy densities, respectively. The corresponding 
Lagrangian density is obtained by changing the sign of 
the potential energy term, and Eq. (1) follows by stan- 
dard variational arguments. This form of the energy 
density leads to two important observations. First, we 
note that the energy density simplifies considerably if 
U describes a soliton and is given by Eq. I jlOj l. Specifi- 
cally, use of the equation of motion allows us to write the 
energy density as £ so \ = u 2 A(u). The specific form of 
Eq. IjlUfl is sufficiently simple that the energy of a soliton 
can be calculated analytically and involves only elemen- 
tary functions. 

It is also useful to consider the total energy associated 
with an arbitrary solution, u(x,t), of Eq. J2J as given 
by the integral over all space of the energy density, £ , 
of Eq. I|ll[) . Recognizing perfect differentials when they 
arise and making use of the equation of motion, Eq. 
we find the expected result that the energy is independent 
of time for an arbitrary choice of u(x, t). (This result 
assumes either that u(x,t) vanishes as \x\ — > oo or that 
it satisfies periodic boundary conditions in x.) It is also 
useful to consider the time dependence of the the integral 
of u over all space, 



U = u(x, t) dx 
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It is clear from the equation of motion that d 2 U/dt 2 can 
be expressed as an integral of perfect differentials. Hence, 
d 2 U/dt 2 — if u vanishes at spatial infinity or is periodic. 
Thus, the time dependence of U is elementary and can 
include only a constant term and a term linear in t. As 
we shall see below, U is independent of time when u(x, t) 
is periodic. 



NUMERICAL CONSIDERATIONS 

We would like to investigate a number of questions 
associated with the stability of the solitons of Eq. (II UP . 



Although the simplicity of the analytic form of these soli- 
tons suggests that it may be possible to solve the problem 
of infinitesimal stability analytically, we have elected to 
consider this problem numerically. To this end, it is con- 
venient to re- write Eq. J3J as two first-order equations. 
We obtain 
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where w = du/dx (and incidentally v = ds/dt). (Note 
that the first of Eqs. I|14l) ensures that the spatial integral 
of u is independent of time if v is chosen to be periodic.) 
Equations (|14f> are well-suited to numerical solution us- 
ing a variant of the two-step Lax-Wendroff method • 
We consider the function, u(x, t), on a primary mesh of 
equally spaced points, (pAx,qAt), where it has the val- 
ues u pq . (Evidently, we must demand Ax < (3 At in order 
to satisfy the usual Courant condition, which is neces- 
sary but not sufficient to ensure numerical stability.) The 
mesh is then extended to include half-integer values of p 
and/or q. It is evident from Eqs. (|14f> and i|15|) that the 
functions u, v, and / live naturally on the points (p, q) 
and (p + 1/2, q + 1/2) and that w lives naturally on the 
points (p + 1/2, q) and (p, q + 1/2). 

The algorithm is implemented as follows: Establish the 
initial values of u and v on the points (p, 0). (If u is a 
solitonic solution, we simply have v(x, 0) = —(3u(x,0).) 
At every time step, the values of u> p +i/2,g are obtained 
as 
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and these values can be used to construct f pq using a sim- 
ilarly symmetric difference formula. One then proceeds 
to the evaluation of the functions u and v a time At/2 
later. In proceeding from q to q + 1/2, we define, 
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The values of 
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Up. q and /p+j./2 g+i/2 are obtained as in 
Eq. ijTljl) . and w p+1 / 2 , g +i/2 is obtained as in Eq. l(T7l) . A 
slightly different procedure is adopted in going from time 
q+ 1/2 to q+ 1. In this case, 

Ax / v 

Up.q+l = U P,q+~^ \ V p+l/2,q+l/2 - «f>-l/2,g+l/2 J ! (18) 

the values of w p +i/2,q+i and f p ,q+i are obtained as in 
Eq. I|16l) . and Vp.g+i is obtained in analogy with Eq. (15). 
This algorithm is both fast and stable in practice. (For 
the periodic boundary conditions and the choice of Ax = 
0.1 and At = 0.001, used below, it was possible to follow 
10 6 time steps without discernible loss of accuracy.) 



4 



It is useful to note that energy of Eq. ljTT|) is not rig- 
orously conserved by this numerical algorithm. In the 
following numerical examples, the energy was found to 
decrease at a roughly constant rate proportional to Ax 2 . 
This fact was used to make an appropriate choice of 
Ax = 0.1. The corresponding value of At = 0.001 was 
selected to yield full numerical stability. 



NUMERICAL RESULTS 
Small amplitude noise 

Our primary numerical concern is to study the stabil- 
ity of the solitonic solutions of Eq. with respect to 
"infinitesimal" perturbations. We employ the parame- 
ters B\ = —16.6 and £>2 = 79.5 adopted in |2J, for which 
flo 0.650. We will show results for an initial soliton 
with velocity (3 = j3\ « 0.735. This soliton has a width 
(i.e., full width at half maximum) of roughly 6.24, which 
is the minimum width possible for the values of B\ and 
i?2 considered. There is, of course, no reason to believe 
that a soliton on a discrete lattice with finite Ax will have 
a profile identical to the analytic form of Eq. i|10[l . The 
use of this analytic form in establishing the initial values 
of u and v thus inevitably introduces a measure of noise 
into the numerical system. Since there is no other "natu- 
ral" choice for the initial form of the solitonic excitation, 
this noise represents the best approximation to infinites- 
imal perturbations that can be realized in a numerical 
study. 

In an analytic approach to the question of infinitesimal 
stability, one considers the time evolution of the sum of 
the soliton under investigation and a small excitation, 
5u(x,t) — ip(x,t). The equation of motion Q is then 
expanded to first order in ijj, and expressed in terms of t 
and £ = x — (3t, 
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It follows that solutions to this (non-Hermitcan) equation 
can be written as the product of functions an d 
exp (Xt). If one or more of the resulting values of A has 
a positive real part, the corresponding ip\(x,t) will grow 
exponentially with time, and the initial solitonic solution 
will be locally unstable. Since it is our aim to detect 
precisely such exponential stabilities (if present), it is of 
no consequence that the numerical noise introduced by 
the finite mesh size is small. Exponential instabilities will 
be apparent if they are present. The finite size of Ax also 
means that there is a smallest wave length perturbation 
which can be studied on the lattice. In practice, potential 
instabilities involving such wave lengths will be invisible 
to numerical studies only if they are orthogonal to those 
wave lengths which can be investigated reliably with the 
Ax chosen. While this is not impossible, it is unlikely. 
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FIG. 2: The difference d(£) between the time-averaged numer- 
ical soliton and the analytic soliton for the minimum width 
soliton with /3 = /3i « 0.735. The average has been performed 
over 1000 units of time, during which the soliton travels more 
than 100 times its own width. 



Results were obtained with Ax = 0.1 and At = 0.001. 
The spatial lattice was chosen to be periodic with length 
100. For (3 = j3\ w 0.735, the exact energy of the soliton 
is 0.0377. The energy of this initial state is smaller by 
1.5 x 10 -6 when calculated on the lattice. Energy is not 
strictly conserved by the numerical algorithm adopted 
but rather decreases linearly with time over the time in- 
tervals considered. In the present case, energy is lost at 
the rate of 7.3 x 10 -9 per unit time. We have followed this 
soliton for times as long as 1000 units, during which the 
soliton can propagate more than 100 times its own width. 
The energy loss is negligible, and there is absolutely no 
indication of instability. (Note that the discrepancy in 
the initial energy is proportional to Ax 2 ; the rate of en- 
ergy loss scales like Ax 3 .) 

We can illustrate soliton stability in the following man- 
ner. We first determine the location of the maximum of 
the soliton as a function of time. The constancy of its 
velocity over large time intervals provides an initial in- 
dication of the stability of the soliton. In the present 
case, this velocity is found to be stable roughly 0.02% 
less than the initial velocity of the analytic soliton. (This 
error scales with Ax 2 .) There are, of course, small fluc- 
tuations in the location of both the maximum density 
and, hence, the velocity due to the presence of noise. For 
the present example, such fluctuations in the location of 
the maximum are never greater than 0.004, which is 25 
times smaller than Ax. (These fluctuations also scale 
like Ax 2 .) Having identified the position of soliton as 
a function of time, each time frame is shifted to locate 
the soliton at a common point. A time-averaged soliton 
is then constructed in order to minimize the effects of 
noise. The difference between the time-averaged soliton 
and the analytic soliton is shown in Figure The peak 
value of the time-averaged soliton is slightly (i.e., roughly 
0.05%) higher than that of the analytic soliton, and it is 



5 



somewhat narrower than its analytic counterpart. (The 
size of these differences again scales with Ax 2 .) This 
demonstrates the claim that the analytic solitons are not 
identical to solitons on a finite mesh. Further, the sys- 
tematic discrepancy between these two solutions is the 
source of and has a magnitude comparable to that of the 
noise in the system. 

We now consider the nature of the "lattice noise" in 
the system as a function of time by subtracting the time 
averaged soliton from the full u(x, t) at each time update 
and constructing the root mean square of the resulting 
noise as a function of time. If the soliton is stable, the 
resulting rms noise should be bounded as a function of 
time. If the soliton were unstable, however, we would ex- 
pect to find systematic differences in the vicinity of the 
soliton maximum which are well above noise level and 
which grow exponentially with time. The spatial dis- 
tribution of noise at later times shows no sign of such 
systematic effects, and its magnitude is the same both 
near and far from the location of the soliton. The cal- 
culated rms noise is shown in Figure 01 as a function of 
time. Again, there is no sign of such instabilities. Since 
qualitatively similar results are found for other values of 
j3, we conclude that the solitons of Eq. © are stable with 
respect to small perturbations. 
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FIG. 3: Time evolution of rms noise level a for the minimal 
width soliton over 1000 units of time. 

It is also possible to study soliton stability in the pres- 
ence of larger amplitude noise. This is most easily done 
by choosing a form of u(x, 0) which consists of both the 
(analytic) soliton of interest and a linear combination 
of the lowest k < K periodic waves on the interval L, 
ak sin (2nkx/L + <pk), with phases chosen at random and 
amplitudes chosen at random subject to a constraint on 
the overall rms noise level at t = 0. The analysis pro- 
ceeds as above. We have considered the case of K = 10 
with an initial rms noise as large as 5% of the maximum 
amplitude of the soliton. The results are similar to those 
found for small amplitude noise: There are no indications 
of soliton instability. 



Soliton genesis 

It is also instructive to consider finite-amplitude distur- 
bances and to see how a localized but non-solitonic ini- 
tial state evolves with time. To illustrate this, we choose 
u(x,0) to be the minimum width soliton of Eq. (|1(J|) . In 
this case, however, we distort the second initial condi- 
tion and choose v(x, 0) = — p(3u(x,0) with p = 0.5. 
Thus, the initial field is not solitonic. The time evo- 
lution shows that this initial pulse "sheds" matter and 
changes its shape through the emission of a smaller soli- 
ton, which moves in the opposite direction, and small 
amplitude waves, which run ahead of the solitons with 
velocity /3 « 1. The two solitons are captured in Fig- 
ure 0| at t = 50. The velocity of the larger soliton is 
(3 = 0.799 and its maximum is at x — 139.515 whereas 
the smaller has (3 = —0.948 and maximum at x = 52.871. 
The shape of each of these solitons is accurately described 
using Eq. l|10fl with the corresponding measured velocity. 
These two solitons account for virtually all of the initial 
energy of the system; approximately 0.3% of this energy 
is associated with the small amplitude motion distinct 
from the solitons. In Figure [5] the two solitons have been 
subtracted out, and only the difference is plotted. This 
confirms that the shapes are indeed solitonic. 
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FIG. 4: A minimal width soliton with an initial velocity, /3, 
50% lower than the corresponding analytic value, shown at 
t — 50. It has divided into two solitons of different sizes, 
propagating in opposite directions. Small- amplitude waves 
run ahead of the solitons with velocity j3 w 1; the region be- 
tween the solitons is essentially noise free. (See also Figure|SJ . 
Note that the length of the periodic lattice has been increased 
here to avoid interference effects between the solitons and the 
leading small amplitude waves. 

Similar results have been obtained for other non- 
solitonic initial pulse forms (e.g., Gaussian pulses). In 
short, for the cases explored, non-solitonic initial exci- 
tations evolve into solitons and small amplitude non- 
solitonic disturbances. In infinite space, dispersion en- 
sures that the solitonic and non-solitonic components will 
become spatially distinct and that the amplitude of the 
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FIG. 5: The graph in Figure 2] with the two solitons sub- 
tracted out to leave only the small-amplitude waves running 
ahead of the solitons. 

latter will decrease with time. This is obviously not the 
case for the periodic lattice considered here. 
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FIG. 6: Decaying soliton (fully drawn) with k = 0.05, initially 
at i = 0. The dashed curves depict the analytic solitons 
with the instantaneous velocity of the numeric solitons. The 
numbers above the peaks indicate the running time, and their 
particular values have been chosen for illustrative purposes. 
The soliton has in fact wrapped around the periodic lattice 
more than 9 times during the time interval of the simulation. 



Solitons and dissipation 

It is also possible to consider the consequences of dissi- 
pation on soliton propagation. The inclusion of viscosity 
in the Navier-Stokes velocity results in an additional term 
on the right of Eq. Q of the form nd 3 u/dx 2 dt. This term 
is readily incorporated in our numerical approach by the 
inclusion of the term +Kdv/d£ in Eq. I)15[l. We have per- 
formed numerical studies with the value k = 0.05. With 
this choice of k, the height of the soliton is reduced by 
roughly 70% at t — 990, has travelled more than 100 
times its initial width. As energy is dissipated, the soli- 
ton accelerates, and its profile changes with the expected 
drop in its amplitude. Over the entire time range con- 
sidered, we find that the soliton profile is consistent with 
the analytic soliton profile of Eq. (|10|) appropriate for the 
corresponding instantaneous velocity of the pulse. This 
is illustrated in Figure 6, which shows the comparison of 
analytic solitons (in infinite space) and these numerical 
results including dissipation at several times. 

For several reasons, this agreement is necessarily only 
approximate. First, some time is required for the soliton 
profile to adjust to the exact form corresponding to its 
instantaneous velocity. Obviously, only a limited time 
is available for this adjustment in the presence of dissi- 
pation. More importantly, the time-independence of the 
spatial integral of u is not affected by the inclusion of dis- 
sipation. Thus, u{x, t) approaches a constant value for all 
x as t — > oo. On a periodic lattice, as here, this constant 
is non-zero. This effect is clearly seen in Figure 6, and it 
is obviously not included in the analytic form of Eq. I|10|) 
valid in infinite space. Figure 6 shows no indication of 
the catastrophic break-up of the soliton into small ampli- 
tude waves which might be anticipated in the presence of 



strong dissipation. It should be noted that magnitude of 
the dissipation considered here is large compared to what 
is to be expected in biomembranes and nerves, where lit- 
tle or no change in pulse shape is observed over distances 
roughly 20 times than the pulse width. 



CONCLUSIONS 

We have considered here a number of tests of the stabil- 
ity of the solitons associated with the modified Boussi- 
nesq equation, Eq. @. After finding the analytic form 
of these solitons, we turned to a numerical investigation 
(with periodic boundary conditions) of the their stabil- 
ity with respect to various perturbations. These solitons 
were found to be stable with respect to the "smallest pos- 
sible" perturbations inevitably induced by the finite size 
of the numerical mesh and to finite but small periodic 
perturbations. Solitons are found to be produced by ar- 
bitrary localized but non-solitonic initial excitations. Fi- 
nally, we have shown that solitons retain their character- 
istic properties even in the presence of relatively strong 
dissipation. It was argued in Ref. that the measured 
compression modulus of lipids of biological membranes 
are suitable for the production of solitons. 

These findings may be of immediate relevance for the 
propagation of the action potential in nerve axons 0. 
The solitons described above are subject to friction and 
dissipation. Nerve membranes are not homogeneous, 
i.e., they vary both in thickness (e.g., at the site of the 
soma) and in the specific composition of lipids and pro- 
teins. Elastic constants may therefore vary locally. In the 
present paper we have shown that neither noise nor dis- 
sipation affect the propagation of solitary waves as such 
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but rather lead only to slight changes in amplitude and 
velocity. These pulses are therefore likely to be robust 
with respect to the unavoidable variance in shape and 
composition of biological membranes and to dissipative 
hydrodynamic processes which accompany the observed 
thickness changes in nerves . Thus, the present results 
suggest that a model of nerve pulses as stable solitons is 
viable even in a realistic physiological environment and 
and that such a model may provide an immediate and re- 
liable explanation of associated mechanical 0] and ther- 
modynamic 0-0-0 effects that remain unexplained in 
the presently accepted Hodgkin-Huxley model 0]. 

We thank Hans Fogedby, Mogens H0gh Jensen, Bo- 
Sture Skagerstam, and Erwin Neher for valuable discus- 
sions. 
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